%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------  Simulate Equilibrium  --------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Setup workspace
close all
clear 

machine = '/Users/giulia/';

set(0,'defaultTextInterpreter','latex')

M = 15;

g_euler = 0.57721;

scale = 10^5;

%% Load matching function estimates

load(path_mf)
options = optimset('Algorithm','interior-point','display','on','MaxIter',90000,'MaxFunEvals',90000,'TolFun',10^-9,'TolX',10^-9,'TolCon',10^-9);

cobb_par = zeros(15,2);

for i =1:15
    col = freights(:,i)~=0;
    s = ships(col,i);
    f = freights(col,i);
    m = matches(col,i);
    lb = [0,0];
    ub = [Inf,0.95];
    theta0 = [rand(1),rand(1)];
    [cobb_par(i,:),~,~] = fmincon(@(theta)cobb(theta,f,s,m),theta0,[],[],[],[],lb,ub,[],options);
end

s_values = [1,linspace(0.95,1.95,10)];
v_values = linspace(0.75,5,11);

clearvars -except s_values v_values scale machine cobb_par

%% Load cost estimates

load(path_ships);

load(path_exports);


%% Efficient allocation

gamma_new = cobb_par(:,2);

[U_eff,m_eff,f_eff,s_eff,J_eff,W_eff,Gdest_eff,Gunc_eff,meanprices_eff,Jf_eff ,ccp_eff ,lambda_eff, lambda_f_eff] = efficient_equil(s_ineff,f_ineff,U_ineff,J_ineff,Jf_ineff,Gdest_ineff,meanprices_ineff,cobb_par,gamma_new,cs,cu,c_inv,sigma,beta,ksi,g_euler,r,delta,K_obs,E,S_obs,ineff,sigma_exp);

w_eff = Welfare(Gdest_eff,Gunc_eff,ccp_eff,s_eff,f_eff,m_eff,cu./sigma,cs./sigma,c_inv,K_obs,ksi,r,E,sigma);
